Cases and social distancing correlation analysis, SCC and SMC

library(tidyverse)
library(plotly)
library(sf)
library(mapview)
library(tigris)
library(censusapi)
library(leaflet)
library(lehdr)
library(usmap)
library(lmtest)
library(pracma)
library(lmtest)
library(forecast)
library(vars)
library(rvest)
library(RSelenium)
library(seleniumPipes)
library(dLagM)

options(
  tigris_class = "sf",
  tigris_use_cache = TRUE
)

Get the case data, first for SCC.

# remDr <- remoteDriver(
#   remoteServerAddr = "192.168.86.25",
#   port = 4445L
# )
# remDr$open()
# 
# remDr$navigate("https://app.powerbigov.us/view?r=eyJrIjoiZTg2MTlhMWQtZWE5OC00ZDI3LWE4NjAtMTU3YWYwZDRlOTNmIiwidCI6IjBhYzMyMDJmLWMzZTktNGY1Ni04MzBkLTAxN2QwOWQxNmIzZiJ9")
# 
# webElem <- remDr$findElements(using = "class", value = "column")
# 
# cases <-
#   1:length(webElem) %>% 
#   map(function(x){
#     webElem[[x]]$getElementAttribute("aria-label") %>% as.character()
#   }) %>% 
#   unlist() %>% 
#   as.data.frame()
# 
# scc_cumulative_cases <-
#   cases %>% 
#   rename(text = ".") %>% 
#   filter(grepl("Total_cases",text)) %>% 
#   separate(text, c("date","cases"), sep = "\\.") %>% 
#   mutate(
#     date = 
#       substr(date,6,nchar(.)) %>% 
#       as.Date("%A, %B %d, %Y"),
#     cases = 
#       substr(cases,13,nchar(.)) %>% 
#       as.numeric()
#   )
# 
# saveRDS(scc_cumulative_cases, file = "/Users/simonespeizer/Documents/2020 Spring Quarter/CEE 218Z/covid19/Simone_Speizer/scc_cumulative_cases.rds")

scc_cumulative_cases <- readRDS("/Users/simonespeizer/Documents/2020 Spring Quarter/CEE 218Z/covid19/Simone_Speizer/scc_cumulative_cases.rds")

Also for SMC.

# remDr$navigate("https://app.powerbigov.us/view?r=eyJrIjoiMWI5NmE5M2ItOTUwMC00NGNmLWEzY2UtOTQyODA1YjQ1NWNlIiwidCI6IjBkZmFmNjM1LWEwNGQtNDhjYy1hN2UzLTZkYTFhZjA4ODNmOSJ9")
# 
# webElem <- remDr$findElements(using = "class", value = "column")
# 
# tests <-
#   1:length(webElem) %>% 
#   map(function(x){
#     webElem[[x]]$getElementAttribute("aria-label") %>% as.character()
#   }) %>% 
#   unlist() %>% 
#   as.data.frame()
# 
# tests_clean <-
#   tests %>% 
#   rename(text = ".") %>% 
#   separate(text, c("date","test_text"), sep = "\\.") %>% 
#   separate(test_text, c(NA,"type",NA,"tests")) %>% 
#   mutate(
#     date = 
#       substr(date,23,nchar(.)) %>% 
#       as.Date("%A, %B %d, %Y"),
#     tests = 
#       tests %>% 
#       as.numeric()
#   ) %>% 
#   spread(
#     key = type,
#     value = tests
#   ) %>% 
#   mutate(
#     total = Negative + Pending + Positive,
#     perc_positive = Positive/total,
#     perc_positive_movavg = movavg(perc_positive, 7, type = "s")
#   )
# 
# smc_cumulative_cases <- tests_clean %>%
#   mutate(cumulative_cases = cumsum(Positive), cumulative_negative = cumsum(Negative), cumulative_total = cumulative_cases+cumulative_negative, perc_positive_cumulative = cumulative_cases*100 / cumulative_total, perc_positive_cumulative_mov = movavg(perc_positive_cumulative, 7, type = "s"))
# 
# saveRDS(smc_cumulative_cases, file = "/Users/simonespeizer/Documents/2020 Spring Quarter/CEE 218Z/covid19/Simone_Speizer/smc_cumulative_cases.rds")

smc_cumulative_cases <- readRDS("/Users/simonespeizer/Documents/2020 Spring Quarter/CEE 218Z/covid19/Simone_Speizer/smc_cumulative_cases.rds")

Get social distancing data.

scc_blockgroups <-
  block_groups("CA","Santa Clara", cb=F, progress_bar=F) %>% 
  st_transform('+proj=longlat +datum=WGS84')

smc_blockgroups <-
  block_groups("CA","San Mateo", cb=F, progress_bar=F) %>% 
  st_transform('+proj=longlat +datum=WGS84')

bay_sd <- readRDS("/Users/simonespeizer/pCloud Drive/Shared/SFBI/Restricted Data Library/Safegraph/covid19analysis/bay_socialdistancing_v2.rds") %>% 
  mutate(date = date_range_start %>%  substr(1,10) %>% as.Date())

# obtaining weekends
weekends <- bay_sd %>% 
  filter(!duplicated(date)) %>% 
  arrange(date) %>% 
  mutate(weekend = ifelse((date %>% as.numeric()) %% 7 %in% c(2,3), T, F)) %>% 
  dplyr::select(date,weekend)

bay_sd <- bay_sd %>% left_join(weekends)

Santa Clara County

scc_cases_sd_daily <- bay_sd %>% 
  filter(origin_census_block_group %in% scc_blockgroups$GEOID) %>%
  group_by(date) %>%
  summarize(total_at_home = sum(completely_home_device_count), total_devices = sum(device_count)) %>%
  mutate(
    percent_at_home = total_at_home*100/total_devices, 
    percent_leaving_home = (100 - percent_at_home),
  ) %>% 
  left_join(
    scc_cumulative_cases
  ) %>% 
  filter(date >= min(scc_cumulative_cases$date))

# get the baseline percent of people at home
pre_case_growth_at_home_scc <- bay_sd %>%
  filter(date < min(scc_cumulative_cases$date)) %>%
  filter(origin_census_block_group %in% scc_blockgroups$GEOID) %>%
  summarize(total_at_home = sum(completely_home_device_count), total_devices = sum(device_count)) %>%
  mutate(percent_at_home = total_at_home*100/total_devices, percent_leaving_home = (100 - percent_at_home))

# include change in percent leaving home
scc_cases_sd_daily <- scc_cases_sd_daily %>%
  mutate(leaving_home_dif = percent_leaving_home - pre_case_growth_at_home_scc$percent_leaving_home[1],
         log_cases = log(cases))

# compute number of differences for stationarity
ndiffs(scc_cases_sd_daily$cases)

## [1] 2

ndiffs(scc_cases_sd_daily$log_cases[-1])

## [1] 2

ndiffs(scc_cases_sd_daily$leaving_home_dif)

## [1] 1

scc_test_big <- scc_cases_sd_daily %>%
  mutate(
    dlog_cases = c(NA,diff(log_cases)),
    d2log_cases = c(NA,NA,diff(log_cases, differences = 2)),
    dcases = c(NA,diff(cases)),
    d2cases = c(NA,NA,diff(dcases, differences = 2)),
    dleaving = c(NA,diff(leaving_home_dif)),
    d2leaving = c(NA,NA,diff(leaving_home_dif, differences = 2)),
    cases_mov = movavg(cases, 7, type = "s"),
    log_cases_mov = movavg(log_cases, 7, type = "s"),
    dlog_cases_mov = c(NA,diff(log_cases_mov)),
    d2log_cases_mov = c(NA,NA,diff(log_cases_mov, differences = 2)),
    dcases_mov = c(NA,diff(cases_mov)),
    d2cases_mov = c(NA,diff(dcases_mov)),
    leaving_mov = movavg(leaving_home_dif, 7, type = "s"),
    dleaving_mov = c(NA,diff(leaving_mov)),
    d2leaving_mov = c(NA,diff(dleaving_mov)),
    leaving4 = c(rep(NA,4), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-4)]),
    dleaving4 = c(NA,diff(leaving4)),
    d2leaving4 = c(NA,NA,diff(leaving4, differences = 2)),
    leaving3 = c(rep(NA,3), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-3)]),
    leaving3_mov = movavg(leaving3, 7, type = "s"),
    dleaving3_mov = c(NA,diff(leaving3_mov)),
    d2leaving3_mov = c(NA,NA,diff(leaving3_mov, differences = 2)),
    leaving4_mov = movavg(leaving4, 7, type = "s"),
    dleaving4_mov = c(NA,diff(leaving4_mov)),
    d2leaving4_mov = c(NA,NA,diff(leaving4_mov, differences = 2)),
    leaving5 = c(rep(NA,5), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-5)]),
    leaving5_mov = movavg(leaving5, 7, type = "s"),
    dleaving5_mov = c(NA,diff(leaving5_mov)),
    d2leaving5_mov = c(NA,NA,diff(leaving5_mov, differences = 2)),
    leaving6 = c(rep(NA,6), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-6)]),
    leaving6_mov = movavg(leaving6, 7, type = "s"),
    dleaving6_mov = c(NA,diff(leaving6_mov)),
    d2leaving6_mov = c(NA,NA,diff(leaving6_mov, differences = 2)),
    leaving7 = c(rep(NA,7), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-7)]),
    leaving7_mov = movavg(leaving7, 7, type = "s"),
    dleaving7_mov = c(NA,diff(leaving7_mov)),
    d2leaving7_mov = c(NA,NA,diff(leaving7_mov, differences = 2)),
    leaving8 = c(rep(NA,8), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-8)]),
    leaving8_mov = movavg(leaving8, 7, type = "s"),
    dleaving8_mov = c(NA,diff(leaving8_mov)),
    d2leaving8_mov = c(NA,NA,diff(leaving8_mov, differences = 2)),
    leaving9 = c(rep(NA,9), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-9)]),
    leaving9_mov = movavg(leaving9, 7, type = "s"),
    dleaving9_mov = c(NA,diff(leaving9_mov)),
    d2leaving9_mov = c(NA,NA,diff(leaving9_mov, differences = 2)),
    leaving10 = c(rep(NA,10), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-10)]),
    leaving10_mov = movavg(leaving10, 7, type = "s"),
    dleaving10_mov = c(NA,diff(leaving10_mov)),
    d2leaving10_mov = c(NA,NA,diff(leaving10_mov, differences = 2)),
    past_cases = c(NA, scc_cases_sd_daily$cases[1:(nrow(scc_cases_sd_daily)-1)]), 
    cases_growth_daily = (dcases / past_cases)*100,
    cases_growth_daily_mov = movavg(cases_growth_daily, 7, type = "s")
  ) %>% 
  filter(date >= "2020-03-01")

ndiffs(scc_cases_sd_daily$cases_growth_daily)

## [1] 0

Tests

Granger

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 3,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:3) + Lags(dleaving_mov, 1:3)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:3)
##   Res.Df Df      F    Pr(>F)    
## 1     60                        
## 2     63 -3 7.7236 0.0001912 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 3,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:3) + Lags(d2log_cases_mov, 1:3)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     60                 
## 2     63 -3 1.4405 0.2399

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 4,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:4) + Lags(dleaving_mov, 1:4)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:4)
##   Res.Df Df      F   Pr(>F)    
## 1     57                       
## 2     61 -4 7.6572 5.21e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 4,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:4) + Lags(d2log_cases_mov, 1:4)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:4)
##   Res.Df Df      F Pr(>F)
## 1     57                 
## 2     61 -4 0.8209 0.5172

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 5,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:5) + Lags(dleaving_mov, 1:5)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:5)
##   Res.Df Df      F    Pr(>F)    
## 1     54                        
## 2     59 -5 5.7937 0.0002365 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 5,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:5) + Lags(d2log_cases_mov, 1:5)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:5)
##   Res.Df Df      F Pr(>F)
## 1     54                 
## 2     59 -5 0.7215 0.6102

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 6,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:6) + Lags(dleaving_mov, 1:6)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:6)
##   Res.Df Df      F    Pr(>F)    
## 1     51                        
## 2     57 -6 4.7436 0.0006565 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 6,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:6) + Lags(d2log_cases_mov, 1:6)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:6)
##   Res.Df Df      F Pr(>F)
## 1     51                 
## 2     57 -6 0.8913 0.5083

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 7,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:7) + Lags(dleaving_mov, 1:7)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:7)
##   Res.Df Df      F  Pr(>F)  
## 1     48                    
## 2     55 -7 2.9577 0.01166 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 7,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:7) + Lags(d2log_cases_mov, 1:7)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:7)
##   Res.Df Df      F Pr(>F)
## 1     48                 
## 2     55 -7 1.6548 0.1431

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 8,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:8) + Lags(dleaving_mov, 1:8)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:8)
##   Res.Df Df      F    Pr(>F)    
## 1     45                        
## 2     53 -8 5.0121 0.0001761 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 8,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:8) + Lags(d2log_cases_mov, 1:8)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:8)
##   Res.Df Df      F  Pr(>F)  
## 1     45                    
## 2     53 -8 1.9809 0.07091 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 9,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:9) + Lags(dleaving_mov, 1:9)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:9)
##   Res.Df Df      F  Pr(>F)   
## 1     42                     
## 2     51 -9 3.8325 0.00133 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 9,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:9) + Lags(d2log_cases_mov, 1:9)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:9)
##   Res.Df Df      F Pr(>F)
## 1     42                 
## 2     51 -9 1.6247 0.1394

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 10,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:10) + Lags(dleaving_mov, 1:10)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:10)
##   Res.Df  Df      F   Pr(>F)   
## 1     39                       
## 2     49 -10 3.6794 0.001578 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 10,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:10) + Lags(d2log_cases_mov, 1:10)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:10)
##   Res.Df  Df     F Pr(>F)
## 1     39                 
## 2     49 -10 1.462 0.1907

Granger for not moving average

grangertest(
  d2log_cases ~ dleaving,
  order = 3,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:3) + Lags(dleaving, 1:3)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:3)
##   Res.Df Df      F   Pr(>F)   
## 1     60                      
## 2     63 -3 4.9045 0.004095 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 3,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:3) + Lags(d2log_cases, 1:3)
## Model 2: dleaving ~ Lags(dleaving, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     60                 
## 2     63 -3 1.3908 0.2543

grangertest(
  d2log_cases ~ dleaving,
  order = 4,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:4) + Lags(dleaving, 1:4)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:4)
##   Res.Df Df      F    Pr(>F)    
## 1     57                        
## 2     61 -4 7.3475 7.628e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 4,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:4) + Lags(d2log_cases, 1:4)
## Model 2: dleaving ~ Lags(dleaving, 1:4)
##   Res.Df Df      F  Pr(>F)  
## 1     57                    
## 2     61 -4 2.3217 0.06763 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  d2log_cases ~ dleaving,
  order = 5,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:5) + Lags(dleaving, 1:5)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:5)
##   Res.Df Df      F    Pr(>F)    
## 1     54                        
## 2     59 -5 5.0124 0.0007638 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 5,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:5) + Lags(d2log_cases, 1:5)
## Model 2: dleaving ~ Lags(dleaving, 1:5)
##   Res.Df Df      F Pr(>F)
## 1     54                 
## 2     59 -5 1.7951 0.1294

grangertest(
  d2log_cases ~ dleaving,
  order = 6,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:6) + Lags(dleaving, 1:6)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:6)
##   Res.Df Df      F  Pr(>F)   
## 1     51                     
## 2     57 -6 3.2324 0.00908 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 6,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:6) + Lags(d2log_cases, 1:6)
## Model 2: dleaving ~ Lags(dleaving, 1:6)
##   Res.Df Df      F Pr(>F)
## 1     51                 
## 2     57 -6 1.3455 0.2546

grangertest(
  d2log_cases ~ dleaving,
  order = 7,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:7) + Lags(dleaving, 1:7)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:7)
##   Res.Df Df      F   Pr(>F)   
## 1     48                      
## 2     55 -7 3.1376 0.008239 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 7,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:7) + Lags(d2log_cases, 1:7)
## Model 2: dleaving ~ Lags(dleaving, 1:7)
##   Res.Df Df      F Pr(>F)
## 1     48                 
## 2     55 -7 1.0694 0.3975

grangertest(
  d2log_cases ~ dleaving,
  order = 8,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:8) + Lags(dleaving, 1:8)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:8)
##   Res.Df Df      F Pr(>F)   
## 1     45                    
## 2     53 -8 3.4128 0.0038 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 8,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:8) + Lags(d2log_cases, 1:8)
## Model 2: dleaving ~ Lags(dleaving, 1:8)
##   Res.Df Df      F Pr(>F)
## 1     45                 
## 2     53 -8 1.6413   0.14

grangertest(
  d2log_cases ~ dleaving,
  order = 9,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:9) + Lags(dleaving, 1:9)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:9)
##   Res.Df Df      F   Pr(>F)   
## 1     42                      
## 2     51 -9 2.9897 0.007637 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 9,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:9) + Lags(d2log_cases, 1:9)
## Model 2: dleaving ~ Lags(dleaving, 1:9)
##   Res.Df Df      F Pr(>F)
## 1     42                 
## 2     51 -9 1.5613 0.1587

grangertest(
  d2log_cases ~ dleaving,
  order = 10,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:10) + Lags(dleaving, 1:10)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:10)
##   Res.Df  Df      F  Pr(>F)  
## 1     39                     
## 2     49 -10 2.3634 0.02696 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 10,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:10) + Lags(d2log_cases, 1:10)
## Model 2: dleaving ~ Lags(dleaving, 1:10)
##   Res.Df  Df      F Pr(>F)
## 1     39                  
## 2     49 -10 0.8308 0.6021

Granger with case growth rate

grangertest(
  cases_growth_daily ~ dleaving,
  order = 2,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:2) + Lags(dleaving, 1:2)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:2)
##   Res.Df Df      F Pr(>F)
## 1     63                 
## 2     65 -2 0.8535 0.4308

grangertest(
  dleaving ~ cases_growth_daily,
  order = 2,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:2) + Lags(cases_growth_daily, 1:2)
## Model 2: dleaving ~ Lags(dleaving, 1:2)
##   Res.Df Df      F  Pr(>F)  
## 1     63                    
## 2     65 -2 3.3287 0.04224 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 3,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:3) + Lags(dleaving, 1:3)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:3)
##   Res.Df Df     F  Pr(>F)  
## 1     60                   
## 2     63 -3 2.343 0.08205 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 3,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:3) + Lags(cases_growth_daily, 1:3)
## Model 2: dleaving ~ Lags(dleaving, 1:3)
##   Res.Df Df      F  Pr(>F)  
## 1     60                    
## 2     63 -3 3.5172 0.02035 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 4,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:4) + Lags(dleaving, 1:4)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:4)
##   Res.Df Df      F   Pr(>F)   
## 1     57                      
## 2     61 -4 4.3608 0.003808 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 4,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:4) + Lags(cases_growth_daily, 1:4)
## Model 2: dleaving ~ Lags(dleaving, 1:4)
##   Res.Df Df      F  Pr(>F)  
## 1     57                    
## 2     61 -4 2.7221 0.03821 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 5,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:5) + Lags(dleaving, 1:5)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:5)
##   Res.Df Df      F   Pr(>F)   
## 1     54                      
## 2     59 -5 3.8974 0.004343 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 5,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:5) + Lags(cases_growth_daily, 1:5)
## Model 2: dleaving ~ Lags(dleaving, 1:5)
##   Res.Df Df      F  Pr(>F)  
## 1     54                    
## 2     59 -5 3.1499 0.01444 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 6,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:6) + Lags(dleaving, 1:6)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:6)
##   Res.Df Df      F    Pr(>F)    
## 1     51                        
## 2     57 -6 4.9432 0.0004699 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 6,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:6) + Lags(cases_growth_daily, 1:6)
## Model 2: dleaving ~ Lags(dleaving, 1:6)
##   Res.Df Df      F  Pr(>F)  
## 1     51                    
## 2     57 -6 2.4259 0.03866 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 7,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:7) + Lags(dleaving, 1:7)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:7)
##   Res.Df Df      F    Pr(>F)    
## 1     48                        
## 2     55 -7 4.3839 0.0007929 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 7,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:7) + Lags(cases_growth_daily, 1:7)
## Model 2: dleaving ~ Lags(dleaving, 1:7)
##   Res.Df Df      F Pr(>F)
## 1     48                 
## 2     55 -7 1.4992 0.1903

grangertest(
  cases_growth_daily ~ dleaving,
  order = 8,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:8) + Lags(dleaving, 1:8)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:8)
##   Res.Df Df      F   Pr(>F)   
## 1     45                      
## 2     53 -8 3.1117 0.006981 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 8,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:8) + Lags(cases_growth_daily, 1:8)
## Model 2: dleaving ~ Lags(dleaving, 1:8)
##   Res.Df Df      F  Pr(>F)  
## 1     45                    
## 2     53 -8 2.9246 0.01022 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 9,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:9) + Lags(dleaving, 1:9)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:9)
##   Res.Df Df    F  Pr(>F)  
## 1     42                  
## 2     51 -9 2.14 0.04708 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 9,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:9) + Lags(cases_growth_daily, 1:9)
## Model 2: dleaving ~ Lags(dleaving, 1:9)
##   Res.Df Df      F  Pr(>F)  
## 1     42                    
## 2     51 -9 2.2421 0.03783 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 10,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:10) + Lags(dleaving, 1:10)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:10)
##   Res.Df  Df      F Pr(>F)  
## 1     39                    
## 2     49 -10 2.0268 0.0567 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 10,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:10) + Lags(cases_growth_daily, 1:10)
## Model 2: dleaving ~ Lags(dleaving, 1:10)
##   Res.Df  Df      F   Pr(>F)   
## 1     39                       
## 2     49 -10 3.3384 0.003229 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Also try ardlDlm

# 4th order lag on x only
testing_ardl4 = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 4, q = 4, remove = list(p = c(0,1,2,3), q=c()))
summary(testing_ardl4)

## 
## Time series regression with "ts" data:
## Start = 5, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0208895 -0.0015317  0.0005453  0.0019138  0.0263403 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0012166  0.0010638  -1.144 0.257296    
## X.4         -0.0007294  0.0013659  -0.534 0.595297    
## Y.1         -0.0499038  0.1301945  -0.383 0.702851    
## Y.2          0.2326294  0.1336737   1.740 0.086937 .  
## Y.3          0.4296713  0.1104830   3.889 0.000255 ***
## Y.4          0.0512616  0.1204743   0.425 0.671996    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.007618 on 60 degrees of freedom
## Multiple R-squared:  0.2415, Adjusted R-squared:  0.1782 
## F-statistic:  3.82 on 5 and 60 DF,  p-value: 0.004534

testing_ardl4 = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 4, q = 4, remove = list(p = c(0,1,2,3), q=c()))
summary(testing_ardl4)

## 
## Time series regression with "ts" data:
## Start = 5, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.107955 -0.014411  0.001501  0.010923  0.095726 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.006614   0.004592  -1.441 0.154912    
## X.4          0.002060   0.001416   1.455 0.150881    
## Y.1         -0.681884   0.126119  -5.407 1.17e-06 ***
## Y.2         -0.571356   0.138753  -4.118 0.000119 ***
## Y.3         -0.441065   0.144300  -3.057 0.003340 ** 
## Y.4         -0.192237   0.105759  -1.818 0.074105 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03626 on 60 degrees of freedom
## Multiple R-squared:  0.3523, Adjusted R-squared:  0.2983 
## F-statistic: 6.526 on 5 and 60 DF,  p-value: 6.593e-05

# 5th order lag on x only
testing_ardl5 = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 5, q = 5, remove = list(p = c(0,1,2,3,4), q=c()))
summary(testing_ardl5)

## 
## Time series regression with "ts" data:
## Start = 6, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.021418 -0.001004  0.001145  0.002117  0.015301 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.0018571  0.0009125  -2.035   0.0464 *
## X.5          0.0012634  0.0011535   1.095   0.2779  
## Y.1          0.1035728  0.1085414   0.954   0.3439  
## Y.2         -0.0214885  0.1097219  -0.196   0.8454  
## Y.3          0.1739958  0.1154550   1.507   0.1372  
## Y.4         -0.0700718  0.1040173  -0.674   0.5032  
## Y.5          0.1172137  0.1014753   1.155   0.2528  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.006404 on 58 degrees of freedom
## Multiple R-squared:  0.2163, Adjusted R-squared:  0.1353 
## F-statistic: 2.669 on 6 and 58 DF,  p-value: 0.0235

testing_ardl5 = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 5, q = 5, remove = list(p = c(0,1,2,3,4), q=c()))
summary(testing_ardl5)

## 
## Time series regression with "ts" data:
## Start = 6, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.102828 -0.010813  0.003215  0.012146  0.081908 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.007695   0.004079  -1.887   0.0642 .  
## X.5         -0.003154   0.001248  -2.527   0.0143 *  
## Y.1         -0.599547   0.111878  -5.359 1.51e-06 ***
## Y.2         -0.303334   0.133295  -2.276   0.0266 *  
## Y.3         -0.319233   0.136188  -2.344   0.0225 *  
## Y.4          0.183904   0.134464   1.368   0.1767    
## Y.5          0.218774   0.094149   2.324   0.0237 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03142 on 58 degrees of freedom
## Multiple R-squared:  0.5057, Adjusted R-squared:  0.4546 
## F-statistic: 9.891 on 6 and 58 DF,  p-value: 1.691e-07

# 6th order lag on x only
testing_ardl6 = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 6, q = 6, remove = list(p = c(0,1,2,3,4,5), q=c()))
summary(testing_ardl6)

## 
## Time series regression with "ts" data:
## Start = 7, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.020128 -0.002107  0.001037  0.001714  0.014682 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.0019052  0.0009642  -1.976   0.0531 .
## X.6          0.0013845  0.0011791   1.174   0.2453  
## Y.1          0.1002644  0.1327362   0.755   0.4532  
## Y.2          0.0306193  0.1104764   0.277   0.7827  
## Y.3          0.2379076  0.1110240   2.143   0.0365 *
## Y.4         -0.1005564  0.1190132  -0.845   0.4018  
## Y.5          0.1042246  0.1055226   0.988   0.3275  
## Y.6         -0.1164147  0.1036735  -1.123   0.2663  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.006468 on 56 degrees of freedom
## Multiple R-squared:  0.2147, Adjusted R-squared:  0.1165 
## F-statistic: 2.187 on 7 and 56 DF,  p-value: 0.04917

testing_ardl6 = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 6, q = 6, remove = list(p = c(0,1,2,3,4,5), q=c()))
summary(testing_ardl6)

## 
## Time series regression with "ts" data:
## Start = 7, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.08285 -0.01097  0.00638  0.01140  0.07022 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0118381  0.0036709  -3.225 0.002105 ** 
## X.6          0.0003069  0.0011401   0.269 0.788801    
## Y.1         -0.6410868  0.1137748  -5.635 5.92e-07 ***
## Y.2         -0.4546921  0.1185650  -3.835 0.000321 ***
## Y.3         -0.5892889  0.1205805  -4.887 8.96e-06 ***
## Y.4         -0.1817903  0.1235876  -1.471 0.146904    
## Y.5         -0.3268229  0.1184621  -2.759 0.007821 ** 
## Y.6         -0.3825815  0.0853473  -4.483 3.70e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02722 on 56 degrees of freedom
## Multiple R-squared:  0.6393, Adjusted R-squared:  0.5942 
## F-statistic: 14.18 on 7 and 56 DF,  p-value: 1.979e-10

# 7th order lag on x only
testing_ardl7 = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 7, q = 7, remove = list(p = c(0,1,2,3,4,5,6), q=c()))
summary(testing_ardl7)

## 
## Time series regression with "ts" data:
## Start = 8, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.022767 -0.001574  0.000301  0.001407  0.011354 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.0009224  0.0008347  -1.105  0.27400   
## X.7          0.0010192  0.0009881   1.031  0.30691   
## Y.1          0.0715185  0.1106610   0.646  0.52083   
## Y.2          0.3304048  0.1104324   2.992  0.00417 **
## Y.3          0.1465600  0.0914905   1.602  0.11501   
## Y.4          0.0871186  0.0955925   0.911  0.36616   
## Y.5          0.1287338  0.0991496   1.298  0.19967   
## Y.6         -0.0561249  0.0880841  -0.637  0.52670   
## Y.7         -0.2286607  0.0867580  -2.636  0.01094 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005351 on 54 degrees of freedom
## Multiple R-squared:  0.3909, Adjusted R-squared:  0.3007 
## F-statistic: 4.332 on 8 and 54 DF,  p-value: 0.0004419

testing_ardl7 = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 7, q = 7, remove = list(p = c(0,1,2,3,4,5,6), q=c()))
summary(testing_ardl7)

## 
## Time series regression with "ts" data:
## Start = 8, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.060278 -0.004423  0.001805  0.007273  0.048478 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.0041398  0.0031885  -1.298  0.19969   
## X.7          0.0020389  0.0009014   2.262  0.02774 * 
## Y.1         -0.3416279  0.1056659  -3.233  0.00209 **
## Y.2         -0.1038290  0.1125298  -0.923  0.36028   
## Y.3         -0.3320478  0.1053011  -3.153  0.00264 **
## Y.4          0.0772674  0.1140214   0.678  0.50088   
## Y.5         -0.0936313  0.0996331  -0.940  0.35152   
## Y.6         -0.1625195  0.0999608  -1.626  0.10981   
## Y.7          0.2141471  0.0787674   2.719  0.00879 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0215 on 54 degrees of freedom
## Multiple R-squared:  0.7712, Adjusted R-squared:  0.7373 
## F-statistic: 22.76 on 8 and 54 DF,  p-value: 9.715e-15

# all
testing_ardl = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 15, q = 15)
summary(testing_ardl)

## 
## Time series regression with "ts" data:
## Start = 16, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -2.986e-03 -9.923e-04  7.876e-05  1.039e-03  3.087e-03 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.0026772  0.0007511  -3.564  0.00165 **
## X.t          0.0016962  0.0011297   1.501  0.14686   
## X.1          0.0006484  0.0013834   0.469  0.64368   
## X.2          0.0013695  0.0012367   1.107  0.27957   
## X.3         -0.0014865  0.0013466  -1.104  0.28105   
## X.4          0.0012411  0.0012553   0.989  0.33309   
## X.5          0.0002892  0.0011818   0.245  0.80884   
## X.6         -0.0002395  0.0011725  -0.204  0.83994   
## X.7          0.0008168  0.0012522   0.652  0.52070   
## X.8         -0.0004827  0.0013803  -0.350  0.72976   
## X.9          0.0004167  0.0011616   0.359  0.72304   
## X.10         0.0033635  0.0011412   2.947  0.00723 **
## X.11        -0.0002858  0.0012217  -0.234  0.81713   
## X.12         0.0004645  0.0012905   0.360  0.72217   
## X.13        -0.0031023  0.0011311  -2.743  0.01159 * 
## X.14         0.0035882  0.0011178   3.210  0.00388 **
## X.15         0.0004624  0.0012469   0.371  0.71413   
## Y.1          0.0356192  0.1567221   0.227  0.82222   
## Y.2          0.0091882  0.1413520   0.065  0.94873   
## Y.3         -0.0156067  0.1392002  -0.112  0.91170   
## Y.4         -0.0622273  0.1420075  -0.438  0.66533   
## Y.5         -0.1140753  0.1147079  -0.994  0.33033   
## Y.6         -0.0138435  0.1120046  -0.124  0.90271   
## Y.7         -0.3005001  0.1192491  -2.520  0.01913 * 
## Y.8         -0.0345115  0.1013978  -0.340  0.73668   
## Y.9         -0.0132726  0.0843804  -0.157  0.87639   
## Y.10        -0.1900936  0.0851933  -2.231  0.03569 * 
## Y.11        -0.0330317  0.0700561  -0.472  0.64172   
## Y.12         0.0175546  0.0664028   0.264  0.79385   
## Y.13         0.0139956  0.0630154   0.222  0.82620   
## Y.14        -0.0201684  0.0639287  -0.315  0.75524   
## Y.15        -0.1572236  0.0655611  -2.398  0.02499 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002179 on 23 degrees of freedom
## Multiple R-squared:  0.9325, Adjusted R-squared:  0.8414 
## F-statistic: 10.24 on 31 and 23 DF,  p-value: 1.206e-07

testing_ardl = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 15, q = 15)
summary(testing_ardl)

## 
## Time series regression with "ts" data:
## Start = 16, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0153804 -0.0046369  0.0002912  0.0038447  0.0200498 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0060194  0.0028369  -2.122 0.044841 *  
## X.t          0.0021201  0.0008117   2.612 0.015588 *  
## X.1          0.0019892  0.0008597   2.314 0.029966 *  
## X.2          0.0035492  0.0009228   3.846 0.000824 ***
## X.3          0.0007635  0.0011123   0.686 0.499316    
## X.4          0.0008054  0.0010050   0.801 0.431115    
## X.5          0.0006726  0.0009732   0.691 0.496397    
## X.6          0.0001779  0.0008291   0.215 0.832030    
## X.7          0.0015223  0.0008115   1.876 0.073423 .  
## X.8          0.0006827  0.0008229   0.830 0.415259    
## X.9          0.0008236  0.0008682   0.949 0.352654    
## X.10         0.0021748  0.0008585   2.533 0.018570 *  
## X.11         0.0017887  0.0008939   2.001 0.057338 .  
## X.12         0.0024315  0.0008839   2.751 0.011385 *  
## X.13         0.0001263  0.0009885   0.128 0.899451    
## X.14         0.0020863  0.0009276   2.249 0.034374 *  
## X.15         0.0006511  0.0009032   0.721 0.478274    
## Y.1         -0.7371592  0.2006844  -3.673 0.001262 ** 
## Y.2         -0.5358158  0.1785719  -3.001 0.006382 ** 
## Y.3         -0.3932115  0.1892215  -2.078 0.049052 *  
## Y.4         -0.2606945  0.1764794  -1.477 0.153186    
## Y.5         -0.2385887  0.1541623  -1.548 0.135358    
## Y.6         -0.0653172  0.1382773  -0.472 0.641119    
## Y.7         -0.1650019  0.1254548  -1.315 0.201395    
## Y.8         -0.2115432  0.1224248  -1.728 0.097404 .  
## Y.9         -0.2025611  0.1114602  -1.817 0.082218 .  
## Y.10        -0.3257234  0.1151342  -2.829 0.009512 ** 
## Y.11        -0.2915754  0.1271660  -2.293 0.031332 *  
## Y.12        -0.1639630  0.1302824  -1.259 0.220829    
## Y.13        -0.0242918  0.1241679  -0.196 0.846613    
## Y.14         0.0824732  0.1067281   0.773 0.447545    
## Y.15         0.0886663  0.0772121   1.148 0.262629    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01069 on 23 degrees of freedom
## Multiple R-squared:  0.8363, Adjusted R-squared:  0.6156 
## F-statistic:  3.79 on 31 and 23 DF,  p-value: 0.0007732

# try with not log of cases
testing_ardl_not_log = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2cases, p = 15, q = 15)
summary(testing_ardl_not_log)

## 
## Time series regression with "ts" data:
## Start = 16, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.8077  -4.6345   0.3506   5.5457  14.4709 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.03807    1.87360   0.020 0.983965    
## X.t          1.06282    0.89591   1.186 0.247616    
## X.1          1.72185    0.84293   2.043 0.052708 .  
## X.2          2.08779    0.76098   2.744 0.011574 *  
## X.3         -0.29649    0.91080  -0.326 0.747724    
## X.4         -1.03824    1.00674  -1.031 0.313130    
## X.5         -1.98815    0.95116  -2.090 0.047845 *  
## X.6         -0.29275    0.84875  -0.345 0.733294    
## X.7          0.77840    0.85529   0.910 0.372208    
## X.8          0.58843    0.81083   0.726 0.475340    
## X.9         -0.37824    0.76945  -0.492 0.627685    
## X.10         0.31816    0.71670   0.444 0.661249    
## X.11         0.27307    0.68208   0.400 0.692594    
## X.12         0.82349    0.66682   1.235 0.229316    
## X.13        -0.74080    0.68498  -1.081 0.290689    
## X.14         2.53881    0.76084   3.337 0.002864 ** 
## X.15         0.21156    0.82757   0.256 0.800505    
## Y.1         -1.47509    0.18745  -7.869 5.68e-08 ***
## Y.2         -1.70247    0.28862  -5.899 5.18e-06 ***
## Y.3         -1.83794    0.36017  -5.103 3.62e-05 ***
## Y.4         -1.94224    0.42495  -4.571 0.000136 ***
## Y.5         -1.73150    0.46087  -3.757 0.001026 ** 
## Y.6         -1.62403    0.46350  -3.504 0.001910 ** 
## Y.7         -1.51701    0.46735  -3.246 0.003564 ** 
## Y.8         -1.47255    0.46274  -3.182 0.004151 ** 
## Y.9         -1.33424    0.45187  -2.953 0.007139 ** 
## Y.10        -1.25388    0.42155  -2.974 0.006785 ** 
## Y.11        -1.04299    0.38103  -2.737 0.011740 *  
## Y.12        -0.85170    0.34327  -2.481 0.020842 *  
## Y.13        -0.58238    0.29787  -1.955 0.062827 .  
## Y.14        -0.21879    0.22868  -0.957 0.348633    
## Y.15        -0.10857    0.14319  -0.758 0.455987    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.23 on 23 degrees of freedom
## Multiple R-squared:  0.9325, Adjusted R-squared:  0.8414 
## F-statistic: 10.24 on 31 and 23 DF,  p-value: 1.207e-07

# pre april 15th
scc_test_pre_april15 <- scc_test_big %>%
  filter(date <=  "2020-04-15")

# all with this dataset
testing_ardl_pre = ardlDlm(x = scc_test_pre_april15$dleaving_mov, y = scc_test_pre_april15$d2log_cases_mov, p = 4, q = 4, remove = list(p = c(0,1,2,3)))
summary(testing_ardl_pre)

## 
## Time series regression with "ts" data:
## Start = 5, End = 46
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0193970 -0.0027041 -0.0002263  0.0023051  0.0263998 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.002201   0.001787  -1.231  0.22617   
## X.4         -0.001325   0.001954  -0.678  0.50215   
## Y.1         -0.053757   0.168624  -0.319  0.75172   
## Y.2          0.248812   0.176541   1.409  0.16731   
## Y.3          0.439239   0.143173   3.068  0.00408 **
## Y.4          0.069857   0.155962   0.448  0.65690   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.009699 on 36 degrees of freedom
## Multiple R-squared:  0.2384, Adjusted R-squared:  0.1327 
## F-statistic: 2.254 on 5 and 36 DF,  p-value: 0.06978

testing_ardl_pre = ardlDlm(x = scc_test_pre_april15$dleaving, y = scc_test_pre_april15$d2log_cases, p = 4, q = 4, remove = list(p = c(0,1,2,3)))
summary(testing_ardl_pre)

## 
## Time series regression with "ts" data:
## Start = 5, End = 46
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.101504 -0.022735 -0.004123  0.018615  0.098219 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.009738   0.007404  -1.315 0.196778    
## X.4          0.002415   0.001969   1.227 0.227944    
## Y.1         -0.697380   0.162307  -4.297 0.000126 ***
## Y.2         -0.597481   0.179501  -3.329 0.002022 ** 
## Y.3         -0.469015   0.187412  -2.503 0.017008 *  
## Y.4         -0.213756   0.138632  -1.542 0.131844    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04614 on 36 degrees of freedom
## Multiple R-squared:  0.3643, Adjusted R-squared:  0.276 
## F-statistic: 4.127 on 5 and 36 DF,  p-value: 0.004587

# ardlBoundOrders
ardlBoundOrders(data = scc_test_big %>% dplyr::select(dleaving, d2log_cases), d2log_cases~dleaving)

## $p
##   dleaving
## 1       14
## 
## $q
## [1] 15
## 
## $Stat.table
##            q = 1     q = 2     q = 3     q = 4     q = 5     q = 6     q = 7
## p = 1  -250.9489 -267.5710 -263.3483 -275.1213 -293.5316 -305.3056 -297.1058
## p = 2  -256.9344 -267.4827 -265.8688 -274.3035 -294.2032 -307.1245 -299.1835
## p = 3  -263.4581 -263.4581 -268.7743 -273.5211 -292.4649 -305.1729 -297.3040
## p = 4  -256.9197 -272.2833 -272.2833 -275.1706 -290.9673 -303.4811 -295.7392
## p = 5  -258.4962 -280.4359 -279.0057 -279.0057 -288.9916 -301.9384 -294.4630
## p = 6  -280.4465 -298.7995 -296.9824 -297.1321 -297.1321 -302.8457 -295.1018
## p = 7  -271.7965 -299.8338 -298.2299 -298.3701 -297.4849 -297.4849 -296.7258
## p = 8  -282.1093 -292.7891 -290.7999 -288.9308 -295.7467 -295.3447 -295.3447
## p = 9  -279.1848 -291.0854 -289.6336 -287.9535 -291.7501 -291.7010 -295.5427
## p = 10 -294.6798 -298.3682 -296.6957 -297.1293 -297.0726 -295.1486 -294.7634
## p = 11 -271.2616 -280.4542 -278.5015 -276.5136 -281.9553 -280.0605 -280.1146
## p = 12 -289.7541 -293.7644 -294.2128 -296.9428 -295.2164 -293.2390 -291.4248
## p = 13 -310.9712 -314.3384 -318.1372 -320.8433 -319.1747 -322.3125 -320.3157
## p = 14 -305.1308 -310.9983 -311.7019 -320.4972 -319.1592 -318.0871 -317.2385
## p = 15 -324.2096 -323.0515 -328.8318 -329.7601 -328.1312 -327.9968 -326.5103
##            q = 8     q = 9    q = 10    q = 11    q = 12    q = 13    q = 14
## p = 1  -304.7931 -297.7748 -293.8927 -290.5643 -295.6733 -293.3990 -298.1992
## p = 2  -306.1365 -299.4384 -293.9848 -290.1470 -294.4886 -292.6269 -296.7515
## p = 3  -305.0713 -298.0177 -293.0737 -288.9128 -294.2433 -292.4169 -295.4953
## p = 4  -303.5975 -297.1456 -291.9681 -287.3530 -292.3349 -290.5310 -294.4090
## p = 5  -304.7066 -299.1953 -292.0180 -288.2497 -291.5839 -289.3112 -294.7581
## p = 6  -303.4262 -298.3603 -291.5864 -287.5411 -291.7405 -290.0525 -297.9255
## p = 7  -302.1107 -297.2501 -289.9498 -285.5484 -290.5526 -291.0481 -298.6609
## p = 8  -303.6832 -297.1590 -290.5779 -284.6231 -288.6746 -289.3063 -296.7745
## p = 9  -295.5427 -297.4785 -291.6063 -286.6418 -288.2629 -292.4720 -297.9880
## p = 10 -293.5849 -293.5849 -293.6999 -292.4071 -291.5141 -293.3995 -301.3351
## p = 11 -283.8958 -285.8242 -285.8242 -296.8939 -301.7199 -303.3599 -318.2587
## p = 12 -290.6325 -289.5466 -300.0136 -300.0136 -299.7320 -304.6793 -318.1959
## p = 13 -318.3452 -320.1624 -321.2812 -329.1462 -329.1462 -327.7065 -325.9072
## p = 14 -315.3132 -317.4421 -318.5133 -324.4961 -324.0695 -324.0695 -325.1359
## p = 15 -325.0579 -330.1588 -329.0686 -333.6106 -340.9263 -341.0082 -341.0082
##           q = 15
## p = 1  -317.9986
## p = 2  -316.0589
## p = 3  -314.5318
## p = 4  -312.8589
## p = 5  -314.6978
## p = 6  -315.5104
## p = 7  -320.9857
## p = 8  -318.9927
## p = 9  -317.3532
## p = 10 -320.8633
## p = 11 -330.8053
## p = 12 -329.2307
## p = 13 -342.8803
## p = 14 -343.2450
## p = 15 -342.2451
## 
## $min.Stat
## [1] -343.245

ARDL test might not require stationarity?

# 4th order lag on x only, original variables
testing_ardl4_2 = ardlDlm(x = scc_test_big$leaving_home_dif, y = scc_test_big$cases, p = 4, q = 4, remove = list(p = c(0,1,2,3), q=c(1,2,3,4)))
summary(testing_ardl4_2)

## 
## Time series regression with "ts" data:
## Start = 5, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.619e-14 -1.373e-15 -3.200e-16  5.980e-16  5.774e-14 
## 
## Coefficients:
##               Estimate Std. Error    t value Pr(>|t|)    
## (Intercept)  6.653e-15  2.920e-15  2.279e+00   0.0261 *  
## X.4         -2.633e-15  2.065e-16 -1.275e+01  < 2e-16 ***
## Y.1          4.761e-16  6.643e-17  7.167e+00 1.11e-09 ***
## Y.0          1.000e+00  6.854e-17  1.459e+16  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.974e-15 on 62 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 2.166e+35 on 3 and 62 DF,  p-value: < 2.2e-16

# ardlBoundOrders
ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, cases), cases~leaving_home_dif)

## $p
##   leaving_home_dif
## 1               15
## 
## $q
## [1] 5
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  551.4315 541.7830 535.8696 530.0526 524.9652 519.4883 505.4713 495.0153
## p = 2  541.2776 542.2030 534.6996 528.1887 523.1593 517.9178 505.7249 495.1092
## p = 3  534.7343 534.7343 536.6990 530.0122 524.8594 519.7008 507.5232 496.6103
## p = 4  528.1732 529.7440 529.7440 530.8711 524.8219 519.0623 507.7080 495.9187
## p = 5  521.6816 522.9820 524.8233 524.8233 526.8219 520.8457 508.9555 495.2984
## p = 6  512.8850 514.6301 516.4811 518.0836 518.0836 519.5213 509.8541 496.7291
## p = 7  506.1375 508.0238 509.6129 511.0533 512.9293 512.9293 511.4407 497.9630
## p = 8  498.6552 500.3424 501.2760 501.6296 502.9350 503.2506 503.2506 495.5341
## p = 9  485.4073 486.0239 486.9155 487.9907 489.8081 491.6142 482.9317 482.9317
## p = 10 476.5659 476.4151 477.7322 477.6431 479.4408 480.9793 476.9776 476.8845
## p = 11 462.1229 462.0401 463.9358 464.6950 466.6923 467.5709 463.7286 463.3197
## p = 12 452.4190 449.4032 451.3701 452.5185 454.4904 456.3255 455.2871 454.5650
## p = 13 446.0101 443.5004 444.7262 441.0506 443.0246 444.7632 446.1498 445.3551
## p = 14 414.3307 413.6259 414.6902 415.8575 414.6262 414.8230 416.5557 417.7452
## p = 15 412.1510 406.2678 406.6905 405.6631 404.3863 405.3680 406.8105 407.5832
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  487.5733 482.6726 474.7761 469.6282 463.6572 455.9885 446.9640
## p = 2  487.4514 481.4059 472.0640 466.7745 460.9846 452.7444 446.1624
## p = 3  489.1777 483.1421 472.1753 466.7051 460.7315 452.3941 446.6445
## p = 4  489.3430 483.7567 473.7310 468.1435 462.0371 452.7884 447.6725
## p = 5  488.0198 482.1261 474.4132 469.3686 463.1432 452.8999 447.4158
## p = 6  489.9574 484.1184 475.9503 470.3785 463.1308 453.2783 448.0971
## p = 7  490.3463 483.2657 474.3034 469.7507 463.7878 454.7004 449.0435
## p = 8  484.6181 475.1708 462.4046 457.7744 452.6782 444.1399 438.7220
## p = 9  484.9306 476.5633 464.3612 459.7471 454.5246 446.0604 440.2099
## p = 10 476.8845 477.5860 466.1511 461.6051 456.0736 446.5651 440.2906
## p = 11 464.3016 464.3016 466.2455 461.6123 456.7817 448.4575 441.8958
## p = 12 456.2021 457.6613 457.6613 459.4294 454.9623 448.3276 439.1265
## p = 13 447.1687 448.1837 450.1686 450.1686 452.0140 443.3848 436.2394
## p = 14 419.6475 420.7660 422.2259 423.2039 423.2039 425.1985 410.9562
## p = 15 409.1554 409.5049 411.4599 411.3893 413.3596 413.3596 409.0664
## 
## $min.Stat
## [1] 404.3863

ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, dcases), dcases~leaving_home_dif)

## $p
##   leaving_home_dif
## 1               15
## 
## $q
## [1] 4
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  563.3203 556.5468 550.6648 543.6184 536.7660 525.5309 512.5341 507.2811
## p = 2  556.9574 558.4840 552.6581 545.6140 538.7639 527.3610 514.4775 509.2262
## p = 3  549.4072 549.4072 551.3550 545.3277 538.0811 526.3433 515.1568 509.6884
## p = 4  543.8399 545.2201 545.2201 546.6075 539.6319 527.2861 516.5984 511.1939
## p = 5  535.6211 536.9059 538.8974 538.8974 539.7245 527.6348 517.8986 512.6875
## p = 6  521.7455 523.4562 525.2409 526.6596 526.6596 523.3985 513.2089 508.4337
## p = 7  511.8901 513.2635 515.1744 516.7433 518.4597 518.4597 514.8977 510.1451
## p = 8  515.9819 516.6476 517.4174 516.9702 516.4896 510.0367 510.0367 511.9807
## p = 9  495.7937 497.7929 498.8826 500.0176 501.4963 495.5237 493.2470 493.2470
## p = 10 486.0202 488.0184 489.8807 490.1027 491.1960 484.7809 485.1142 485.6264
## p = 11 476.8587 477.9075 479.9013 480.1035 481.7970 476.3355 475.3342 475.8643
## p = 12 466.2053 466.9407 467.8379 469.6922 471.1108 467.6380 467.1936 468.2445
## p = 13 448.7376 450.7376 452.4398 450.8206 452.6976 451.3775 452.9695 452.4121
## p = 14 430.0209 431.8791 433.0050 423.6393 425.2960 426.3064 428.3032 428.3932
## p = 15 423.2949 425.1817 427.0081 421.2607 423.2438 421.3601 423.3568 421.9097
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  500.9525 495.3731 488.5888 483.6054 478.3645 470.1307 462.2284
## p = 2  502.8943 497.3640 490.5714 485.6030 480.3645 472.0792 464.2245
## p = 3  503.0829 497.3656 491.5017 486.5515 481.0187 473.0953 464.8591
## p = 4  503.9822 496.9930 490.0851 484.9604 479.5173 472.2849 463.1397
## p = 5  505.6865 498.1247 489.4702 484.4219 478.6026 470.9423 463.6693
## p = 6  502.7776 495.0681 483.5585 476.0396 469.9647 461.4657 454.9348
## p = 7  504.6644 497.0398 485.5517 477.8507 471.0057 463.2649 456.8309
## p = 8  506.3176 498.7452 486.2678 478.4237 472.3820 461.7054 454.7676
## p = 9  495.1670 488.6104 476.0440 469.5490 461.5124 454.3281 446.5764
## p = 10 485.6264 485.9962 474.3965 466.8136 457.1185 450.2643 442.2470
## p = 11 477.8021 477.8021 476.3958 468.8136 458.5311 451.7515 444.2388
## p = 12 470.2362 471.3915 471.3915 470.2410 460.2490 453.2161 444.4992
## p = 13 454.4119 454.3142 455.4411 455.4411 451.4153 444.5477 436.7992
## p = 14 430.3769 430.0190 431.9277 432.4159 432.4159 432.6255 422.5280
## p = 15 423.8915 423.3048 425.1451 423.5634 423.5652 423.5652 424.4842
## 
## $min.Stat
## [1] 421.2607

ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, d2cases), d2cases~leaving_home_dif)

## $p
##   leaving_home_dif
## 1               14
## 
## $q
## [1] 15
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  596.8029 585.1344 574.8355 567.6410 534.8172 527.7619 520.8212 515.8249
## p = 2  590.1137 587.0190 576.8176 569.6169 536.7224 529.5562 522.6747 517.6850
## p = 3  576.8215 576.8215 575.6664 568.6428 537.9695 530.9054 524.3234 519.3164
## p = 4  575.9275 569.5017 569.5017 570.2982 539.7050 532.7799 526.2488 521.2463
## p = 5  570.1679 566.1116 564.1554 564.1554 541.6558 534.6894 528.2101 523.2170
## p = 6  558.7235 554.8029 554.6769 534.3003 534.3003 534.5655 527.6211 522.5500
## p = 7  546.1128 542.9943 540.9415 539.9070 528.8434 528.8434 529.0404 524.0259
## p = 8  548.2126 543.6979 540.9000 541.3655 522.3490 522.3181 522.3181 524.2777
## p = 9  535.6439 530.0268 526.6979 528.2749 510.9973 512.6614 513.7313 513.7313
## p = 10 529.4588 524.6606 521.2560 522.7175 506.5407 508.4563 508.6884 510.4387
## p = 11 517.4364 513.5564 508.7137 510.1027 493.8151 495.5684 496.0371 498.0363
## p = 12 506.4488 502.5069 499.8967 501.0570 484.9590 486.6526 485.8424 487.8103
## p = 13 481.8881 480.3362 479.9163 478.9870 471.8291 473.7712 473.8086 475.6352
## p = 14 463.0507 463.3476 454.9262 450.7764 448.0149 448.6843 445.7275 447.4288
## p = 15 445.4927 446.3352 438.4082 440.4080 437.1953 438.5188 436.6688 438.3851
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  509.7835 504.5044 498.9399 486.8007 479.8890 473.8559 467.4932
## p = 2  511.6966 506.4000 500.7271 488.6346 481.2843 475.3213 469.2805
## p = 3  513.1026 507.7074 502.3867 490.1641 482.5489 476.9465 470.6609
## p = 4  514.8333 509.3236 503.7799 491.3077 483.1360 476.8671 471.0562
## p = 5  516.8324 511.3215 505.7229 492.9721 484.5546 477.8199 472.5858
## p = 6  516.8858 511.6751 505.8550 492.9150 484.2152 476.3484 471.2303
## p = 7  518.0338 512.7477 506.2194 492.8402 484.1382 475.5451 470.3127
## p = 8  518.0536 512.3625 504.3034 488.8966 478.0287 464.8852 460.1567
## p = 9  515.3571 509.4286 502.4835 487.2160 475.5037 463.7814 458.4289
## p = 10 510.4387 511.3092 504.2173 488.8262 477.4583 465.1542 459.7079
## p = 11 499.4619 499.4619 500.6853 485.8143 475.4783 464.9744 458.4922
## p = 12 489.5813 490.9129 490.9129 480.0151 469.4743 460.5481 452.5700
## p = 13 477.3555 478.8883 470.2526 470.2526 469.8552 460.5796 452.6569
## p = 14 446.0042 446.6970 444.8195 439.8721 439.8721 439.9391 430.7741
## p = 15 436.0841 437.8856 438.1497 431.7956 432.0592 432.0592 432.7404
## 
## $min.Stat
## [1] 430.7741

ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, d2log_cases), d2log_cases~leaving_home_dif)

## $p
##   leaving_home_dif
## 1               15
## 
## $q
## [1] 15
## 
## $Stat.table
##            q = 1     q = 2     q = 3     q = 4     q = 5     q = 6     q = 7
## p = 1  -247.1158 -262.6566 -259.7280 -271.1864 -289.5665 -306.0241 -297.9282
## p = 2  -260.1151 -270.6964 -271.7295 -276.1488 -293.3896 -305.4478 -297.7216
## p = 3  -272.1176 -272.1176 -275.2289 -276.2051 -295.3421 -308.4182 -300.1463
## p = 4  -257.3371 -276.7596 -276.7596 -275.8190 -293.9966 -306.7898 -298.4602
## p = 5  -257.5284 -277.3622 -275.5188 -275.5188 -292.2041 -304.8192 -296.4861
## p = 6  -273.3595 -291.6139 -290.3658 -291.6970 -291.6970 -303.0178 -294.7003
## p = 7  -272.3023 -300.7645 -299.2209 -297.8259 -297.7265 -297.7265 -296.0480
## p = 8  -282.0520 -292.6506 -290.6616 -288.8317 -296.8763 -295.5956 -295.5956
## p = 9  -278.5346 -293.5507 -292.1824 -290.5130 -294.2584 -292.7547 -295.4100
## p = 10 -292.4685 -297.2686 -295.9552 -295.7238 -295.5189 -293.5769 -292.4145
## p = 11 -271.5680 -281.0769 -279.1167 -277.1916 -282.6667 -280.7441 -280.4900
## p = 12 -297.7511 -299.3264 -302.0489 -303.4335 -301.4789 -299.6070 -297.8830
## p = 13 -294.5457 -294.4147 -296.4900 -301.1529 -299.8009 -299.7930 -299.7014
## p = 14 -306.9080 -314.9017 -318.8187 -328.3627 -326.4117 -325.5119 -323.5907
## p = 15 -322.5124 -320.8542 -327.1413 -330.1459 -329.1888 -328.1993 -326.7715
##            q = 8     q = 9    q = 10    q = 11    q = 12    q = 13    q = 14
## p = 1  -303.1134 -296.7856 -293.6299 -289.3689 -299.7575 -300.6864 -295.6658
## p = 2  -302.9945 -295.8307 -292.4216 -288.5985 -301.3512 -307.0429 -306.3532
## p = 3  -304.6415 -297.5186 -292.3243 -288.1629 -299.7060 -306.6731 -307.4808
## p = 4  -303.5515 -296.1163 -291.4783 -286.9199 -299.7326 -306.3824 -308.1563
## p = 5  -302.3817 -295.2922 -290.2753 -285.3606 -297.7328 -304.5140 -306.6937
## p = 6  -302.8127 -297.2628 -290.6908 -286.4765 -298.3092 -305.8850 -308.5185
## p = 7  -301.7097 -296.3767 -289.7787 -285.6182 -297.0538 -304.7369 -307.9607
## p = 8  -300.4574 -295.3230 -288.0325 -283.6450 -295.5137 -305.0929 -309.1502
## p = 9  -295.4100 -295.2668 -288.7073 -282.8341 -294.0422 -303.4526 -307.4371
## p = 10 -291.3980 -291.3980 -289.6450 -284.6432 -292.1702 -303.2680 -306.9536
## p = 11 -283.4290 -285.4861 -285.4861 -290.4073 -294.3568 -302.0437 -309.2631
## p = 12 -297.7796 -295.7838 -302.8100 -302.8100 -305.2507 -312.4671 -322.4358
## p = 13 -298.3913 -296.9694 -298.4937 -302.1586 -302.1586 -314.0261 -322.6509
## p = 14 -322.4463 -325.2010 -324.5372 -329.7972 -328.1825 -328.1825 -331.6931
## p = 15 -325.3757 -330.3256 -329.2409 -333.4078 -340.3607 -341.1173 -341.1173
##           q = 15
## p = 1  -311.9700
## p = 2  -317.2822
## p = 3  -315.7170
## p = 4  -313.7580
## p = 5  -311.9534
## p = 6  -314.2362
## p = 7  -314.2484
## p = 8  -319.3696
## p = 9  -317.3698
## p = 10 -315.6900
## p = 11 -319.0116
## p = 12 -328.8293
## p = 13 -327.2857
## p = 14 -340.9265
## p = 15 -341.4775
## 
## $min.Stat
## [1] -341.4775

ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, cases_growth_daily), cases_growth_daily~leaving_home_dif)

## $p
##   leaving_home_dif
## 1               15
## 
## $q
## [1] 15
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  393.1757 375.3039 357.0862 340.7902 334.0467 282.5013 275.6036 269.8120
## p = 2  379.1824 376.6348 352.4453 339.6024 331.6322 284.2639 277.5686 271.7322
## p = 3  369.7678 369.7678 353.5871 340.5742 332.7598 285.2172 277.9381 272.7965
## p = 4  341.3875 343.1355 343.1355 340.2579 334.4932 286.7561 279.8589 274.3915
## p = 5  347.9692 347.9300 334.7346 334.7346 329.4404 288.4827 281.6367 276.3697
## p = 6  331.4535 332.4995 313.9106 312.9811 312.9811 290.0421 282.5161 276.0033
## p = 7  296.1998 295.4014 282.5097 284.3730 284.3179 284.3179 284.0541 277.7899
## p = 8  289.1624 290.9319 285.0433 286.7694 288.5234 276.8694 276.8694 278.8337
## p = 9  289.7759 291.0523 276.4356 278.4331 279.7860 271.4070 273.2216 273.2216
## p = 10 276.5552 278.0002 266.3790 267.5950 267.8485 257.1157 259.0252 260.6031
## p = 11 262.4497 263.5248 261.6312 262.3240 264.3200 261.7908 263.7794 264.2309
## p = 12 221.5408 221.4510 222.8320 223.6527 221.4774 223.4773 223.8759 225.7048
## p = 13 209.1545 206.8361 207.7154 208.5987 202.8317 204.6966 201.2008 198.7951
## p = 14 198.2758 198.5349 189.7918 185.3501 175.6170 177.6005 179.4515 180.1959
## p = 15 191.6988 185.1561 186.7437 180.1345 175.2661 175.9340 176.3101 177.5464
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  264.4046 256.5817 251.1818 218.0678 203.8873 198.1708 182.5250
## p = 2  266.4033 258.4472 253.0561 219.7598 204.2247 197.6611 180.6552
## p = 3  266.7027 259.5050 254.8955 221.6009 206.1549 198.7904 181.0040
## p = 4  268.4987 261.2078 256.5957 223.2449 208.1431 200.2166 182.9999
## p = 5  270.1325 262.8564 258.5765 224.8884 209.9436 201.5634 184.0590
## p = 6  269.1356 260.2955 256.8728 220.7890 199.4582 192.4838 175.7583
## p = 7  270.9138 261.6299 258.8486 222.7881 200.2446 190.4380 175.9913
## p = 8  272.7483 263.6298 260.8486 224.7564 202.1607 192.4373 176.3502
## p = 9  267.5247 259.9905 258.4628 225.1660 200.5352 189.2311 177.3842
## p = 10 260.6031 261.2215 258.7955 226.5752 201.8393 190.2866 179.1897
## p = 11 260.3329 260.3329 259.4055 227.6676 203.3513 191.5684 180.2022
## p = 12 224.9154 226.8497 226.8497 225.0855 194.6972 188.4177 171.0937
## p = 13 197.6100 199.4856 199.3437 199.3437 196.3416 188.1092 172.8223
## p = 14 170.8064 171.9388 173.2853 168.8885 168.8885 170.6939 160.8284
## p = 15 179.4827 177.8003 175.8217 170.5812 159.8029 159.8029 161.4884
## 
## $min.Stat
## [1] 159.8029

testing_ardl4_growth = ardlDlm(x = scc_test_big$leaving_home_dif, y = scc_test_big$cases_growth_daily, p = 5, q = 5, remove = list(p = c(0,1,2,3)))
summary(testing_ardl4_growth)

## 
## Time series regression with "ts" data:
## Start = 6, End = 70
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.7931  -1.1739  -0.0644   1.1413  12.4465 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.77637    4.04403   0.687 0.495158    
## X.4          0.28624    0.15813   1.810 0.075547 .  
## X.5         -0.19359    0.16047  -1.206 0.232641    
## Y.1          0.23420    0.11525   2.032 0.046815 *  
## Y.2          0.19142    0.11495   1.665 0.101355    
## Y.3         -0.09892    0.11362  -0.871 0.387578    
## Y.4          0.37583    0.09750   3.855 0.000297 ***
## Y.5          0.08068    0.09769   0.826 0.412309    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.755 on 57 degrees of freedom
## Multiple R-squared:  0.8358, Adjusted R-squared:  0.8156 
## F-statistic: 41.43 on 7 and 57 DF,  p-value: < 2.2e-16

Scatter plot to look at thresholds

scc_nolag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving_mov) %>%
  cbind(`Lag` = "0")

scc_3lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving3_mov) %>%
  cbind(`Lag` = "3")
colnames(scc_3lag)[ncol(scc_3lag)-1] = "leaving_mov"

scc_4lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving4_mov) %>%
  cbind(`Lag` = "4")
colnames(scc_4lag)[ncol(scc_4lag)-1] = "leaving_mov"

scc_5lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving5_mov) %>%
  cbind(`Lag` = "5")
colnames(scc_5lag)[ncol(scc_5lag)-1] = "leaving_mov"

scc_6lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving6_mov) %>%
  cbind(`Lag` = "6")
colnames(scc_6lag)[ncol(scc_6lag)-1] = "leaving_mov"

scc_7lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving7_mov) %>%
  cbind(`Lag` = "7")
colnames(scc_7lag)[ncol(scc_7lag)-1] = "leaving_mov"

scc_8lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving8_mov) %>%
  cbind(`Lag` = "8")
colnames(scc_8lag)[ncol(scc_8lag)-1] = "leaving_mov"

scc_9lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving9_mov) %>%
  cbind(`Lag` = "9")
colnames(scc_9lag)[ncol(scc_9lag)-1] = "leaving_mov"

scc_10lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving10_mov) %>%
  cbind(`Lag` = "10")
colnames(scc_10lag)[ncol(scc_10lag)-1] = "leaving_mov"

scc_scatter <- rbind(scc_nolag, scc_3lag, scc_4lag, scc_5lag, scc_6lag, scc_7lag, scc_8lag, scc_9lag, scc_10lag)

fig_cases_growth <- 
  plot_ly (scc_scatter) %>%
    add_trace(
      x = ~leaving_mov, 
      y = ~cases_growth_daily_mov, 
      frame = ~`Lag`, 
      type = 'scatter', 
      mode = 'markers', 
      showlegend = F
    ) %>%
  animation_button(visible = F) %>%
  layout(xaxis = list(title = '% difference in leaving home'), yaxis = list(title = 'daily case growth rate'), title = list(title = "SCC"))

fig_cases_growth

fig_cases_abs <- 
  plot_ly (scc_scatter) %>%
    add_trace(
      x = ~leaving_mov, 
      y = ~dcases,
      frame = ~`Lag`, 
      type = 'scatter', 
      mode = 'markers', 
      showlegend = F
    ) %>%
  animation_button(visible = F) %>%
  layout(xaxis = list(title = '% difference in leaving home'), yaxis = list(title = 'daily new cases'), title = list(title = "SCC"))

fig_cases_abs

Plot the case growth rate and social distancing data

4 and 6 day lags seem the best, depending on which Granger test metric was used.

scc_test_big %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving4_mov, color="Leaving home")) +
  geom_line(aes(y = cases_growth_daily_mov-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%), 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 4 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "Santa Clara County Case Growth and Social Distancing, 4 Day Lag")

scc_test_big %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving4_mov, color="Leaving home")) +
  geom_line(aes(y = dcases_mov-40, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+40, name = "Daily new cases, 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 4 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "Santa Clara County Case Growth and Social Distancing, 4 Day Lag")

scc_test_big %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving4_mov, color="Leaving home")) +
  geom_line(aes(y = d2cases_mov-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Change in change in cases, 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 4 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "Santa Clara County Case Growth and Social Distancing, 4 Day Lag")

# no moving average
scc_test_big %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving4, color="Leaving home")) +
  geom_line(aes(y = cases_growth_daily-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%)")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 4 days ago relative to before", x = "Date", color = "Data", title = "Santa Clara County Case Growth and Social Distancing, 4 Day Lag")

# 6 day
scc_test_big %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6_mov, color="Leaving home")) +
  geom_line(aes(y = cases_growth_daily_mov-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%), 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "Santa Clara County Case Growth and Social Distancing, 6 Day Lag")

# no moving average
scc_test_big %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6, color="Leaving home")) +
  geom_line(aes(y = cases_growth_daily-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%)")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before", x = "Date", color = "Data", title = "Santa Clara County Case Growth and Social Distancing, 6 Day Lag")

Try to animate.

second_axis <- list(overlaying = "y",
  side = "right",
  title = "Daily case growth rate (%), 7 day moving average")

fig_animate_lines_scc <- 
  plot_ly (scc_scatter %>% mutate(date_numeric = as.numeric(date))) %>%
  add_trace(
    x = ~date_numeric,
    y = ~leaving_mov,
    frame = ~`Lag`,
    type = 'scatter',
    mode = 'marker'
    #name = "% leaving"
  ) %>%
  animation_button(visible = F) %>%
  add_trace(
      x = ~date_numeric, 
      y = ~cases_growth_daily_mov,
      frame = ~`Lag`,
      type = 'scatter', 
      mode = 'lines', 
      yaxis = "y2"
      #name = "cases"
    ) %>%
  layout(xaxis = list(title = 'date, numeric'), yaxis = list(title = 'Change in % of devices leaving home relative to before, 7 day moving average'), yaxis2 = second_axis, title = list(title = "SCC"))

fig_animate_lines_scc

San Mateo County

For now, just plotting the metric of cumulative percent positive

smc_cases_sd_daily <- bay_sd %>% 
  filter(origin_census_block_group %in% smc_blockgroups$GEOID) %>%
  group_by(date) %>%
  summarize(total_at_home = sum(completely_home_device_count), total_devices = sum(device_count)) %>%
  mutate(
    percent_at_home = total_at_home*100/total_devices, 
    percent_leaving_home = (100 - percent_at_home),
  ) %>% 
  left_join(
    smc_cumulative_cases
  ) %>% 
  filter(date >= (min(smc_cumulative_cases$date)-10))

# get the baseline percent of people at home
pre_case_growth_at_home_smc <- bay_sd %>%
  filter(date < min(smc_cumulative_cases$date)) %>%
  filter(origin_census_block_group %in% smc_blockgroups$GEOID) %>%
  summarize(total_at_home = sum(completely_home_device_count), total_devices = sum(device_count)) %>%
  mutate(percent_at_home = total_at_home*100/total_devices, percent_leaving_home = (100 - percent_at_home))

# include change in percent leaving home
smc_cases_sd_daily <- smc_cases_sd_daily %>%
  mutate(leaving_home_dif = percent_leaving_home - pre_case_growth_at_home_smc$percent_leaving_home[1],
         log_cases = log(cumulative_cases))

smc_test_small <- smc_cases_sd_daily %>%
  mutate(
    leaving_mov = movavg(leaving_home_dif, 7, type = "s"),
    leaving4 = c(rep(NA,4), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-4)]),
    leaving4_mov = movavg(leaving4, 7, type = "s"),
    leaving5 = c(rep(NA,5), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-5)]),
    leaving5_mov = movavg(leaving5, 7, type = "s"),
    leaving6 = c(rep(NA,6), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-6)]),
    leaving6_mov = movavg(leaving6, 7, type = "s"),
    leaving7 = c(rep(NA,7), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-7)]),
    leaving7_mov = movavg(leaving7, 7, type = "s"),
    leaving8 = c(rep(NA,8), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-8)]),
    leaving8_mov = movavg(leaving8, 7, type = "s"),
    leaving9 = c(rep(NA,9), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-9)]),
    leaving9_mov = movavg(leaving9, 7, type = "s"),
    dcases = c(NA,diff(cumulative_cases)),
    leaving10 = c(rep(NA,10), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-10)]),
    leaving10_mov = movavg(leaving10, 7, type = "s"),
    dcases = c(NA,diff(cumulative_cases)),
    past_cases = c(NA, smc_cases_sd_daily$cumulative_cases[1:(nrow(smc_cases_sd_daily)-1)]), 
    cases_growth_daily = (dcases / past_cases)*100,
    cases_growth_daily_mov = movavg(cases_growth_daily, 7, type = "s")) %>% 
  filter(date >= "2020-03-01")

smc_test_small %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6, color="Leaving home")) +
  geom_line(aes(y = cases_growth_daily-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%)")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before", x = "Date", color = "Data", title = "San Mateo County Case Growth and Social Distancing, 6 Day Lag")

# moving average
smc_test_small %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6_mov, color="Leaving home")) +
  geom_line(aes(y = cases_growth_daily_mov-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%), 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "San Mateo County Case Growth and Social Distancing, 6 Day Lag")

# try with percent positive cumulatively
smc_test_small %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6, color="Leaving home")) +
  geom_line(aes(y = perc_positive_cumulative-40, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+40, name = "Cumulative percent of tests that are positive")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before", x = "Date", color = "Data", title = "San Mateo County Positive Tests and Social Distancing, 6 Day Lag")

smc_test_small %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6_mov, color="Leaving home")) +
  geom_line(aes(y = perc_positive_cumulative_mov-40, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+40, name = "Cumulative percent of tests that are positive, 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "San Mateo County Positive Tests and Social Distancing, 6 Day Lag")

smc_test_small %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6_mov, color="Leaving home")) +
  geom_line(aes(y = perc_positive_movavg*100-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1-30/100, name = "Percent of tests that are positive, 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "San Mateo County Positive Tests and Social Distancing, 6 Day Lag")

smc_nolag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving_mov) %>%
  cbind(`Lag` = "0")

smc_4lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving4_mov) %>%
  cbind(`Lag` = "4")
colnames(smc_4lag)[2] = "leaving_mov"

smc_5lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving5_mov) %>%
  cbind(`Lag` = "5")
colnames(smc_5lag)[2] = "leaving_mov"

smc_6lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving6_mov) %>%
  cbind(`Lag` = "6")
colnames(smc_6lag)[2] = "leaving_mov"

smc_7lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving7_mov) %>%
  cbind(`Lag` = "7")
colnames(smc_7lag)[2] = "leaving_mov"

smc_8lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving8_mov) %>%
  cbind(`Lag` = "8")
colnames(smc_8lag)[2] = "leaving_mov"

smc_9lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving9_mov) %>%
  cbind(`Lag` = "9")
colnames(smc_9lag)[2] = "leaving_mov"

smc_10lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving10_mov) %>%
  cbind(`Lag` = "10")
colnames(smc_10lag)[2] = "leaving_mov"

smc_scatter <- rbind(smc_nolag, smc_4lag, smc_5lag, smc_6lag, smc_7lag, smc_8lag, smc_9lag, smc_10lag)

fig_cases_growth_smc <- 
  plot_ly (smc_scatter) %>%
    add_trace(
      x = ~leaving_mov, 
      y = ~cases_growth_daily_mov, 
      frame = ~`Lag`, 
      type = 'scatter', 
      mode = 'markers', 
      showlegend = F
    ) %>%
  animation_button(visible = F) %>%
  layout(xaxis = list(title = '% difference in leaving home'), yaxis = list(title = 'daily case growth rate'), title = list(title = "SMC"))

fig_cases_growth_smc